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Abstract 

We introduce a method and algorithm for computing the weighted Moore- 
Penrose inverse of multiple-variable polynomial matrix and the related algo- 
rithm which is appropriated for sparse polynomial matrices. These methods 
and algorithms are generalizations of algorithms developed in |24] to multiple 
variable rational and polynomial matrices and improvements of these algo- 
rithms on sparse matrices. Also, these methods are generalizations of the 
partitioning method for computing the Moore-Penrose inverse of rational and 
polynomial matrices introduced in [25] and |23| to the case of weighted Moore- 
Penrose inverse. Algorithms are implemented in the symbolic computational 
package MATHEMATICA. 

AMS Subj. Class.: 15A09, 68Q40. 

Key words: Weighted Moore-Penrose inverse; rational matrices, polynomial 
matrices; sparse matrices; symbolic computation. 

1 Introduction 

Let C mx ™ be the set of m x n complex matrices, and C™ x ™ is the set of m x n 
complex matrices of rank r: C™ x " = {X G C" IX ™ | rank(A) = r}. For any matrix 
A E c mxn and positive definite Hermitian matrices M and N of the order m and n 
respectively, consider the following equations in X, where * denotes conjugate and 
transpose: 

(1) AX A = A (2) XAX = X 

(3) (MAX)* = MAX (4) {NX A)* = NX A. 

The matrix X satisfying these equations is called the weighted Moore-Penrose in- 
verse of A, and it is denoted by X — A' MN . In the partial case M — I m , N = I n , 
the matrix X = A\ [N comes to the Moore-Penrose inverse A^ of A. 
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As usual, C[si, ■ ■ ■ ,s p ] (resp. C(si, . . . , s p )) denotes the polynomials (resp. ra- 
tional functions) with complex coefficients in the variables si, . . . , s p . The matrices 
of format mxn with elements in C[si, . . . , s p ] (resp. C(si, . . . , s p )) are denoted 
by C[si, . . . , s p ] mxn (resp C(si, . . . , s p ) mxn ). By 7 it is denoted an appropriate 
identity matrix. 

Computation of the Moore-Penrose inverse of one variable polynomial and/or 
rational matrices, based on the Leverrier-Faddeev algorithm, is investigated in 
Q]i[2]>[3j,[I],[I],[6]. Implementation of this algorithm in the symbolic computational 
language MAPLE, is described in [3]. Algorithm for computing the Moore-Penrose in- 
verse of two- variable rational and polynomial matrix is introduced in [7 . A quicker 
and less memory-expensive Effective algorithm for computing the Moore-Penrose 
inverse of one- variable and two- variable polynomial matrix, with respect to those in- 
troduced in [4] and [7] , is presented in [8] . This algorithm is efficient when elements 
of the input matrix are polynomials with only few nonzero addends. 

Papers [9] .[10] .[5] deal with a computation of the Drazin inverse. A general- 
ization of these algorithms, introduced in [11], generates the wide class of outer 
inverses of a rational or polynomial matrix. 

An interpolation algorithm for computing the Moore-Penrose inverse of a given 
one-variable polynomial matrix, based on the Leverrier-Faddeev method, is pre- 
sented in [12] . Algorithms for computing the Moore-Penrose and the Drazin inverse 
of one- variable polynomial matrices based on the evaluation-interpolation technique 
and the Fast Fourier transform are introduced in [13] . Corresponding algorithms 
for two- variable polynomial matrices are introduced in [14] . 

In this paper we consider the set of rational and polynomial matrices and various 
variants of the partitioning method for computing generalized inverses. Grevile's 
partitioning method for numerical computation of generalized inverses is introduced 
in [15]. Two different proofs for Greville's method were presented in [16], [17]. A 
simple derivation of the Grevile's result has been given by Udwadia and Kalaba 
[T5] , In [H] Fan and Kalaba used the approach of determination of the Moore- 
Penrose inverse of matrices using dynamic programming and Belman's principle 
of optimality. Wang in [20] generalizes Grevile's method to the weighted Moore- 
Penrose inverse. 

Many numerical algorithms for computing the Moore-Penrose inverse lack nu- 
merical stability. The Greville's algorithm requires more operations and conse- 
quently it accumulates more rounding errors (see for example [H]). Moreover, it is 
well-known that the Moore-Penrose inverse is not necessarily a continuous function 
of the elements of the matrix. The existence of this discontinuity present further 
problems in the pseudoinverse computation. It is therefore clear that cumulative 
round off errors should totally eliminated, which is possible only by means of the 
symbolic implementation. In the symbolic implementation variables are stored in 
the "exact" form or can be left "unassigned" (without numerical values), resulting 
in no loss of accuracy during the calculation [4]. 

An algorithm for computing the Moore-Penrose inverse of one-variable poly- 
nomial and/or rational matrices, based on the Grevile's partitioning algorithm, is 
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introduced in [22J . An extension of results from [22J to the set of two- variable ratio- 
nal and polynomial matrices is introduced in the paper [23] . In our recent paper |24j 
we propose an algorithm for computing the weighted Moore-Penrose of one- variable 
rational and polynomial matrix. In this work we generalized the results from |24j 
in the following two ways: 

- extends algorithms from [24] to the set of multi-variable rational and polyno- 
mial matrices with complex coefficients, 

- make algorithms from [24) more effective on sparse matrices with a relatively 
small number of nonzero elements. 

The structure of the paper is as follows. In the second section we extend the algo- 
rithm for computing the weighted Moore-Penrose from [20] to the set of multiple- 
variable rational matrices with complex coefficients. Main results are given in the 
third and the fourth section. In Section 3 we adapt previous algorithm to the set 
of polynomial matrices. In the fourth section we consider two effective structures 
which exploit only nonzero addends in polynomial matrices and improve previous 
results on the set of sparse matrices. In the last section we presented an illustrative 
example and compared various algorithms. 

2 Weighted Moore-Penrose inverse for multi- 
variable rational matrices 

Let A(si, . . . , s p ) be complex rational matrix. For the sake of simplicity, we will 
introduce new variables S2 P +i-i = si- Also we will denote the vector of all variables 
si, . . . , S2 P by S = (s\, . . . , S2 P ) and further we will denote A(s\, . . . , s p ) as A(S). 

By Ai(S) we denote the submatrix of A(S) consisting of its first i columns, and 
by a,i(S) is denoted the i-th column of A(S); 

A i (S) = [A i _ 1 (S)\a i (S)],i = 2,...,n, A^S) = ai (S) (2.1) 

We will consider positive definite Hermitian matrices M(S) € C(S) mxm and N(S) S 
C(S) nxn . The leading principal submatrix N t (S) £ C(5) iXl of N(S) is partitioned 

as 

Ni~x(S) h(S) 
l*(S) nu(S) 



Ni(S) 



, i = 2,...,n, (2.2) 



where k(S) £ C(5) ( * and nu(S) is the complex polynomial. By N 1 (S) we 

denote the polynomial rin (S). 

In the following lemma we generalize the representations of the weighted Moore- 
Penrose inverse from |17j.|24) to the set of rational matrices of multiple complex 
variables C(5)" iX ™. 

For the sake of simplicity, by Xi(S) we denote the weighted Moore-Penrose 
inverse corresponding to M(S) and submatrices Ai(S), Ni(S): Xi(S) — Ai(Sy MN , 
for each i = 2, . . . , n. Similarly X±(S) = ai(S) M N 
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Lemma 2.1. Let A(S)eC(S) m * n , assume that M(S)eC(S) mxm , N(S)eC(S) nyn 
are positive definite Hermitian matrices, and let Ai{S) be the submatrix of A(S) 
consisting of its first i columns, as it is defined in <\2. Jp . Assume that the leading 
principal submatrix Ni(S) £ C(S) txl is partitioned as in \2. 2\ . Then the matrices 
Xi(S) can be computed in this way: 

x ren _ / (a* 1 (S)M(S)a 1 (S)y 1 a* 1 (S)M(S), oi(S) + 0, , , 



X^ 1 (S)-(d i (S) + (I-X i _ 1 (S)A^ 1 (S))Nr\(S)l i (S))b^(S) 

b*dS) 

i = 2,...,n, 



where the vectors di(S), Ci(S) and b*(S) are defined by 
di(S) = Xi-^aiiS) 

d(S) = a l (S)-A l _ 1 (S)d l (S) = (I-A l _ 1 (S)X l _ 1 (S))a l (S) 



(c*(S)M(S) Cl (S)) 1 c*(S)M(S), c t (S)^0 

5^ 1 (S)(d* i (S)N i _ 1 (S)-k(S)*)X^ 1 (S), c i (S) = 0, 
and where in b*(S) is 



(2.4) 

(2.5) 
(2.6) 

(2.7) 



Si(S) = n u (S) + d*(S)N i . 1 (S)d i (S) ~ (d*(S)k(S) + l*(S)di(S)) 
-l*(S) (7-X i _ 1 (5)A i _ 1 (5))iVr 1 1 (5)^(S). 



(2.8) 



Also in [T7] authors used a block representation of the inverse Nj 1 (S), which 
we also generalized to the set of rational matrices. 

Lemma 2.2. Let Ni(S) be the partitioned matrix defined in $2.2$ . Assume that 
Ni(S) and Ni^i(S) are both nonsingular. Then 



where 



2Vi_i(S) k(S) 
l*(S) n u (S) 



~] -l 



E^S) MS) 
f*(S) hu(S) 



i = 2, 



,n, 



(2.9) 



h u (S) = {n u (S)-l*(S)Nr_\( S )k(S)y 
MS) = -h»(S)Nr\(S)k(S) 
Ei^(S) = Nr_\(S) + h^(S)MS)f*(S). 



(2.10) 

(2.H) 
(2.12) 
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In view of Lemma 12.11 and Lemma 12.21 respectively, we present the following 
algorithms for computing the weighted Moore-Penrose inverse and the inverse ma- 
trix N^(S) G C(S) IX1 . These algorithms are generalizations of corresponding 
algorithms from [24j to the set of multiple- variable rational matrices with complex 
coefficients. 

Algorithm 2.1. Input: A(S) G C(S) mxn and positive definite matrices M(S) 6 
C(S) mxm and N(S)eC{S) nxn . 

Step 1. Initial value: Compute Xi{S) — a\{Sy defined in \2,ty . 

Step 2. Recursive step: For each i = 2, . . . , n compute Xi(S) performing the 
following four steps: 

Step 2.1. Compute di(S) using <\2.5\ . 

Step 2.2. Compute Ci(S) using \2. 6'|) . 

Step 2.3. Compute b*{S) by means of §2J§ and (TQ|) . 

Step 2.4. Applying \2.4\ compute Xi(S). 

Step 3. The stopping criterion: i — n. Return X n (S). 



Algorithm 2.2. Let JVj(5) = 



be the leading principal subma- 



N^x(S) k(S) 
l*(S) n tl (S) 

trix of positive definite matrix N £ C(S) nxn . Then the inverse matrix N~ 1 (S) can 
be computed as follows: 

Step 1. Initial values: A^ 1 " 1 (S') = n^(S). 

Step 2. Recursive step: For i = 2, . . . , n perform the following steps: 

Step 2.1. Compute hu(S) using <\2.10\i . 

Step 2.2. Compute fi(S) using <\2.11\ . 

Step 2.3. Compute E i ^ 1 (S) using (2A2} . 

Step 2.4. Compute Nf 1 ^) using (pO|) . 
Step 3. For i — n return the inverse matrix N^ 1 (S) = iV^ 1 (>!?). 

We used MATHEMATICA function Together in order to enable simplifications of 
rational expressions (this function joins rational addends together and cancels com- 
mon multipliers in numerator and denominator). 



3 Weighted Moore-Penrose inverse for multi- 
variable polynomial matrices 

Now suppose that A(S) 6 C[5] mx " is multi- variable polynomial matrix. We can 
represent it in the following polynomial form: 

di d 2p Q 

A(S) = ]T A iu ..., i2p si 1 ---s^=j2A I S I , (3.1) 

i 1= i 2p =0 1=0 
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where I=(ii, . . . , i2p), Ai=Ai 1> ... j i 2p are constant mxn matrices, S 1 = s| 1 s 1 ^ ■ • ■ s^fi 
Q = (di, . . . , d,2p) — degA(S). Here di is the degree of the matrix polynomial with 
respect to the variable Sj in A(S). 

If by J we denote J = (j2p, ■ ■ ■ > ji), where J — (ji, . . . > J2p) then it can be easily 

Q 

checked that holds A* (5) = £ A}S J . 

An application of Algorithm 2.1 to the multiple- variable polynomial matrix A(S) 
gives the following result. 

Theorem 3.1. Let us consider A(S) G C[S , ]™ IX " of the form jp and positive 
definite Hermitian matrices M{S) GC(S) mxm and A(S) £ C(5) nx ™. Assume i/iai 
i/ie leading principal submatrix Ni(S) G C(5) 4Xl of N(S) is partitioned as in $2.2$ . 
Then the weighted Moore-Penrose inverse A* MN \S) S C* xm [S l ] corresponding to 
the first i columns in A(S) is of the form 

X i (S)=Al Ni (S) = ^, i = l,...,n, (3.2) 

where Zi(S) G C mXi [S] andYi(S) G C[S], can be computed from Zi-i(S), Yi_i(S), 
Ai^i(S) and a,i(S) using exact recurrence relations. 

Proof. We will prove theorem by the induction. In the case i — 1 exact relations 
for Zx(S) and Y\(S) can be derived from ([2~3]) : 

ax(S) = Ai(S) = ^ Zi(S) = 0, Yx(S) = l 

ai (S) = Ai(5) ^ = a*(S)M(S), Yi(5) = a* 1 (S)M(S)a 1 (S) 

Consider now the inductive step. From the inductive hypothesis we can write 

z ( s) 

Xi-i(S) = y' Hs) ■ Then Xi(S) can be computed by using Step 2 of algorithm 
2.1. From steps 2.1 and 2.2 we have: 

a (q\ y fen (q\ Z t _ 1 (S)a l (S) Di(S) 



Ci(S) = Oi(5) - A i ^ 1 (S)d i {S) = 



Yi-i(S) Yi-i(S) 

Oi(5)yi_i(5) - Ai-iWDiiS) _ d(S) 



Yi^(S) Yi-tiS) 
If Ci(S) ^ 0, according to the Step 2.3 of Algorithm 2.1 we have: 

i[ } cpj^MwgUfc C*(S)M(S)C i (S) Wi(S) 
Otherwise, we need first to evaluate the expression 8i(S). From (|2.8I) we obtain: 
= Tiii(5) + Ni-i{S)-, 



Yf^iS) ^ 'Y^(S) 
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Here we used the inductive hypothesis together with temporary polynomial matrix 
cf>i(S) S C[S'](^ 1 ) xl and polynomial ^(5) are defined by: 



(I - Xi^fflAi-^S)) N-_\(S)h(S) 



Yi-x^Ni-^SMS) - Z^ 1 (S)A t - 1 (S)N^ 1 (S)h(S) MS) 



(3.4) 



y,_i(s)jv,_i(s) 



MS)' 



Also, we use Nf_\(S) = where NdS) £ CfS]**- 1 )***- 1 ) and N^S) e C[S] 

are defined in the next theorem. By collecting addends under the same denominator 
in (|3.3j) we can write Si(S) in the form: 



S t (S) 



A,; (S) 



MS) 

where: 

Ai(S) = n u (S)N l MS)Y;_ 1 (S)Y t - 1 (S)+N l - 1 (S)D*(S)N l MS)D l (S) 

-(K i _ 1 (5)c*(5)/ 4 (5)+r;_ 1 (5)A(^)/*(5))iv l - 1 (5)-/*(^)^(^)^-i(^) 
Ai(s) = yjij (5)ri_i (s)JVi-i (5). 

Now we apply Step 2.3 in the case Ci(S) = and evaluate b*(S): 

Zi-i(S) 



a 4 (5) vi7-i(sy 



Vi-i(S) 



aU(s) (DtWN^iS) - ycmsVKS)) z^s) _ v t (s) 



\(S) 

Let us rewrite now expression (|2.4p in following way: 

/ Di(s) ms) \ Vi(s) ■ 
\Yt-ks) + Ms)) mis) 

Vi(S) 



Wi(S) 



Xi(S) 



Zi-i(S) 



1 



Wi(S) 

Wi 



(S)fti-i(S)Zi-i(S)-(Di{S)Ni-i(S) + <lH(S)) Vi(S) 
N t MS)Y t - 1 (S)V l (S) 



Wi(s)Ms) 

[From the last expression we obviously have that holds: 

_ W i (S)iVi_ 1 (5)^_ 1 (5) - (A(S)^i-i(5) + MS)) Vi(S) 



Zi 



N i _ 1 {S)Y i - 1 (S)V i (S) 









_MS)_ 



Yi = W t (S)MS)- 
This completes the proof of the theorem. □ 
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Theorem 3.2. Let the leading principal submatrix Ni(S) of the positive definite 
matrix N(S) S C[s] nx ™ is partitioned as in $2.2$ . Then the inverse N~ 1 (S) is of 
the form: 



Ni(S) 



F*(S) H U (S) 



Ni(S) Ni(S) 

where £*_i(5) G C^- 1 )*^ 1 ), F*(S) G C^ 1 )* 1 and scalar H U (S) G C[s] can be 
computed from AT i _ 1 (,S l ), l(S), na{S), Ni-±(S) and Ni-i(S) using exact recurrence 
relations. 

Proof. As in the proof of the previous theorem we will use induction and lemma 
2.2 (algorithm 2.2). The case i = 1 is again trivial and we have: 

Ni(S) = l, N 1 (S) = n n (S) 

Let us consider now the inductive step and suppose that N^\(S) = j^' -1 ^ . From 
the relation (|2.10|) we have: 



1 - n u (S)-lUS) N J- 1 [ S ) li(S) (3.5) 



Ni-t^nuiS) - Z*(S)AVi(S)Z 4 (S) Hi(S) 



represent fi(S) in following way 



Therefore, we can write Hu(S) = g , ' . Using the relation (|2.11|) we can 



Nj-^S) Nj^jS) _ ll^N^jS) = Fj(S) 

MS) fti-i(S) Hi(S) MS)' 

Furthermore using the fact that Ni-i(S) is symmetric and positive definite, we 
can conclude that F*(S) = Ni-i(S)k(S) which further implies that: 

f* (S ) = gigj = ^-i(S)k(S) 
H*(S) Hi(S) 

We also used that Hi(S) = H*(S) which can be easily proven from (13.51) . From 
(|2. 12|) we can conclude: 

E = N^jS) H t (S) F t (S)F*(S) 
1 1 Ni-x(S) Ni-i(S) Hi(S) Hi(S) 
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Finally, we can represent 1 (S) in the following matrix form: 

Ei-i(S) 



iV 4 _i(S)Hi(Sf) 

Hi(S) 
1 



Nj-i(S) 
H«(S) 

Ei-!(S) Ni_ x (S)F i (S) 



Nj(S) 

Ms)' 



This completes proof of the theorem. □ 

Now it is easy to construct corresponding algorithms from the theorems 3.1 and 

3.2. 



4 Effective method 

In practice we often work with polynomial matrices A(S) with a relatively small 
number of nonzero coefficients. In that case, previous algorithm is not effective 
because of many operations are redundant. To avoid this problem we will construct 
two appropriate sparse structures for the representation of the polynomial matrix 
A(S) and corresponding effective algorithm for computing A' MN (S). The first sparse 
representation is denoted by Eff and its improvement by Eff', while the second 
structure is denoted by Ef . 

The main idea in the first considered sparse structure is to exploit only non-zero 
coefficient matrices Aj = Ai lt ___^ 2 ^ of the polynomial matrix A(S) given in the 
form ([33]) . 

Definition 4.1. The effective sparse structure of the polynomial matrix A(S), de- 
fined in <\3.1\ . is equal to: 

(4.1) 



Eff A = {(J,Aj)\Aj^O, 0< J <deg A(S)}. 

Also define the index set of this effective structure by: 

lnd A = {J | Aj + 0, < J < deg A(S)} . 

Define operations +, — , • and * on sparse structures as: 

Eff A + Eff B = ES A+B , Eff^ - Eff B = Eff 
Eff A ■ Eff B = Eff A B , Eff^ = Eff^ . 

Denote by eA = \ESa\ = |Ind^| the size of the structure Eff^. 
Obviously we have 

A(S) ■ B(S) = A i B J sI+J > 

I g lnd A 
J € Ind s 



(4.2) 



(4.3) 
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where 

qi+J _ Ji+h j2p+j2 P 
d — h x ■ ■ ■ b 2p 

If C(S) = A(S)B{S) then the elements of Eff c are pairs (K, C K ) where C K is 
defined as the following sum of matrix products: 

C K = a i b k-i (4.4) 

/ e IndA, 
K -I e Ind B 

where Ck ^ 0. Therefore holds ec < &a + &b and Effc = Eff a • Effs can be 
computed in the time 0(ca ■ es). 

Similarly holds for computing the sum C(S) = A(S) + B(S). Elements of Eff c 
are pairs (K, Ck) where values C'k are defined by 

f A K , A K ^0,B K = 

C K ={ B K , B K ^0,A K = (4.5) 

[ Ak + Bk, Ak^Q.B k ^Q 

and satisfy Ck ^ 0. As in the previous case we can conclude that ec < maxjeA, &b} 
and Eff q can be computed in time C^maxje^, e^}). 

Index sets corresponding to addition and multiplication of sparse matrices are 
equal to: 

IndA+s = Ind^i U Ind^, Ind^B = Ind,4 + Inds 
In view of gSJ), we compute Eff^={(J, A}) | (J,Aj)eEff A } in time 0(e A )- 

Usually, coefficient matrices A\ in the polynomial representation (|3.1[) . i.e. in 
the sparse representation (|4.1j) are sparse. Using this fact we can significantly im- 
prove our sparse structure Eff by using an appropriate structure for these constant 
coefficient matrices. 

Definition 4.2. For the constant matrix A = [dij] G C mx ", denote the following 
sparse structure: 

Sp A = {(i,j,a ij )\a ij ^0} (4.6) 
Denote by s A = | Sp^ | the size of the structure Sp A . 

Similarly as in the case of Eff^, we can define elementary operations on these 
sparse structures: 

Sp A + Sp B = {(i,j,aij + bij) | (i,j,a,ij) £ Sp^ V (i,j,bij) G Sp B , aij+bij 
s Pa ■ s Pb = {(*> 3> (kj) | Cjj ajkhj 0, (i, k,a ik ) eSp A A (k, j, b kj ) eSp B } 
Sp^ = {0' ? i » a *j) I (iJ>aij) 6 s Pa} 
In this way, we have the following improvement of the structure Eff : 

Eff^ = {(J, Sp Aj ) \Aj^0,0<J< deg A(S)} (4.7) 
= {(J,{h3,(Aj)uKAjh ;^Q}) \AjjtQ, 0< J < deg A(S)}. 
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It can be seen that the complexity of computing Sp^ + Sp B is 0(sa + s B ) and for 
Sp^ is O(sa)- In the case of multiplication the complexity depends on concrete 
implementation. Suppose that A G c mxn and B £ C nxp . If the triples are sorted 
lexicographically in Sp^ and in Sp B then for every (i,k, Aik) € Sp^ we need to 
find all (k,j,Bkj) € Sp B , i.e. all triples in Sp B which begin by k. If we denote this 
number by s B : 

= \{(k,j,b kj )eSp B \j = i,..., P }\ 

then the complexity of multiplication Sp^ • Sp B is: 

ol sP+m-p). (4.8) 

\(i,k,a ik )eSp A J 

The last addend in (4.8) comes from the fact that we need to construct the sparse 
structure Sp c for the matrix C = AB 6 C mx P. 

We implemented the sparse structure Sp in MATHEMATICA as the structure 
SparseArray. Mathematica offers a sparse representation for matrices, vectors, 
and tensors with SparseArray [25], |26j . Both of the expressions 

SparseArray[{ii, ji}-> Vi, {i 2 , h}~> ^2, ■ •-,] 

SparseArray[{{ii, ji}, {i 2 , j 2 }, ■ ■ •}-> {vi,v 2 , ■ ■ ■}} 
represent the sparse array with elements in positions {ik,jk} having values Vk- 
Operations on sparse matrices are all equivalent to the operations on dense matrices 
[25] . |26j : Plus(+) for matrix addition, Dot (. ) for matrix multiplication, Times (*) 
for multiplication by scalar, etc. 

Therefore, in our implementation we have 

Eff^ = {(J, SparseArray[Aj]) | Aj ^ 0, < J < deg A(S)} . 

Shown fact that basic operations are the same for dense and sparse matrices allows 
us to use the same procedures for basic operations on Eff in cases when Sp is 
embedded in Eff and when it is not. In procedural programming languages we can 
decide to use Sp or not in the beginning of algorithm, depending of the structure 
of input matrices ^4(5*), M(S) and N(S). Similarly, it is possible to change the 
choice of one between these two variants of the structure Eff during the algorithm 
implementation. 

In the second type of the sparse structure for polynomial matrices we represent 
the matrix A(S) in the form A(S) = [aij(S)}, where ciij(S) are scalar polynomials, 
and construct effective sparse structures Eff aij for each aij(S). Effective structure 
Eff a for the scalar polynomial a(S) — ai gi j s defmed similarly as in the 

matrix case (|4.ip : 



Eff„ = {(J,aj) | aj f 0, 1 < J < dega(S)}. 
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Such sparse representation we denote by Ef, and have Ef^ = [Eff .-]. If we use 
notations ef^ = S?=i e oijj then the complexity for the addition is 



E E e «« + e ^ = + effi )- 

\i=l 3=1 / 

After the notations row(_B, k) = 53?=i e fcfcj an d col(A, fc) = J^ i=1 e 0ifc we conclude 
that the complexity of the matrix multiplication C = AB is equal to 

(m n p \ / n m \ 

E E E e «^ e ^ = E E e ^ row (A fc ) 
i=l fc=l j"=l / \fe=l i=l / 



0[ ^ col(A, fc)row(B, fc) 



\fc=i / 

for the multiplication. 

Polynomials in MATHEMATICA are represented in the internal form using the little 
modified Ef sparse structure. For example, two- variable polynomial p(s±,S2) — 
4sf S2 + s 2 + s i s 2 + 3fiiS2 + s 2 + 2s| + 3si + 10 is represented in the following 
MATHEMATICA internal form: 

Plus[ 
10, 

Times [3 , si] , 

Times [2, Power [si, 2]], 

s2, 

Times [3, Power [si, 3], s2] , 

Times [Power [si, 2], Power [s2, 2]], 

Power [s2, 3] , 

Times[4, Power [si, 9], Power[s2, 10]] 
] • 

The last expression is obtained by using MATHEMATICA function FullForm [E] which 
returns an internal representation of the expression E [25) , [26] . This internal form 
of the polynomial p(S), at the top level is the list with length ef p with the head 
Plus. Each element of this list contains the exponent J = (j'1,,7'2) and the value 
p.j (values ji — 0, 1 and j'2 = 0, 1 and are not shown), hence the length of each 
element is 0(1). Also the size of whole structure is 0(ef p ( s )). Therefore, we can 
use this natural polynomial representation in MATHEMATICA and built-in elementary 
operators to implement the effective partitioning method using Ef structure. The 
complexity of these built-in operations are the same as corresponding operations 
defined for Ef structure. 

The next algorithm is the effective partitioning method for computing the we- 
ighted Moore-Penrose inverse of polynomial matrices, suitable for sparse matrices. 
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Generally, the same method can be used with both two presented sparse structures. 
Therefore, we will denote general sparse structure with £ , which can be exchanged 
either by Eff or Ef . Also by O we will denote the general effective structure of an 
appropriate zero matrix. We will use the same symbol for the effective structure of 
the number 0. 

Algorithm 4.1. [Computing the weighted Moore-Penrose inverse A(S) M , S j jy(s) 
of sparse matrix A(S)). 

Input: Effective structures of matrices A(S), M(S), N(S). 

Step 1. In the case £ ai ^ O compute initial values: 



If £ ai = O, then set £z t — O and £y l — £\, where £\ is the corresponding 
sparse structure of the number 1 . 

Step 2. Recursive step: For i = 2, . . . , n perform the following steps 

Step 2.1 Compute: £d i = £z % - x • £a t 

Step 2.2 Compute: £ Ci — £ ai ■ £y i - 1 — • £di 

Step 2.3 If £ Ci 7^ O then compute £y t and £y/ i using 



t Z 1 - o ai ■ t, M , o Yl - c ai ■ t M ■ o, 




Otherwise use the following formulae: 




where the structures Aj and Aj are defined by: 




£d» ■ £ii — £i* ■ £di) + £d' ■ £n^—\ ■ 



£d,) ■ 



We used sparse representations for temporary variables ipi and tpi, defined 




Step 2.4- Now compute £z t and £y t using: 



£y t — £i/> t ■ £\Vi 
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Structures £@ t and £^ i are defined by: 



£&i — £z t -! ■ ^N i _ 1 • — £<ii ■ ^ N i l • £v, - £<ei£vi 

£q/ i = £ii) t ■ £vi 



If we use Ef or EfF sparse structure, £z t is equal respectively to: 

T?f - l" Efei 
Efz * " [Ef* 4 _ 

(j, (eo j ) e EfF 0i , (i, j ) g Kfr„, ; 
(j,(e 4 ),) eEff ei ,(*i)j ■ =o| 



EfF 



z t : 



J, 



o 



(4.9) 



U 





mi 



(9,), =0,(i, (*,),) G EfF*, 



S'iep 2.5. Find f/ie polynomials Zi(S) and Yi(S) from its effective structures and 
compute: 

*m = fff (4.io) 

Cancel the common multipliers in numerator Zi(S) and denominator 
Yi(S), recompute (if necessary) effective structures and continue with 
the next i. 

Step 3. The stopping criterion is i = n. In this case is A' M ^ N ^(S) = X n (S). 

Similarly we can derive a modification of the method introduced in Theorem 3.2 
for computing the inverse matrix N~ 1 (S) in the polynomial form: 



Nj(S) 

Ms) 



(4.11) 



Algorithm 4.2. (Effective computation of N i 1 (S), for i = 1, . . . ,n). 

Input: Effective structure of positive definite Hermitian polynomial matrix N(S) 
of the order n. Notations are the same as in Theorem 3.2. 

Step 1. Generate initial values: N\ = I and N\ = mi and corresponding effective 
structures. 

Step 2. Recursive step: For i = 2, . . . , n perform following steps: 
Step 2.1. Compute: £f [ , — £ nu ■ £ft._ 1 — £* f ■ £ft. 1 ■ £\ i . 
Step 2.2. Compute: £p, — —£jj. • £u- 
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Step 2.3. Compute: £g = £^ 
Step 2.4-. Generate: 

Ni{S) = 
Ni{S) = 



- £f 4 • ■ 

Ei-i(S) Ni-i{S) ■ Fi(S) 

N^^-F^S) Nl^S) 

Ni-i(S) ■ Hi(S) 



As in the previous algorithm, we have also two different representations 
for EF and EfF sparse structures. These relations are similar to OTP| . 



Step 3. Stop criterion for i = n. Inverse matrix N k 1 (S), for every k — 1, . . . ,n is 
equal to: 



N k (S) 



(4.12) 



5 Examples 



We implemented algorithms 2.1, 2.2, 4.1 and 4.2 in the programming language 
MATHEMATICA. An implementation of the EfF sparse structure is also made. Func- 
tions WPolyEf and WPolyEf f implement Algorithm 4.1 using respectively Ef and 
EfF sparse strucure. All basic operations for EfF sparse structure (functions Add, 
Sub, Muls, Mul and TE corresponding to the addition, subtraction, multiplication by 
scalar, multiplication and conjugate-transposion respectively) are also implemented. 

Example 5.1. Let us find the weighted Moore-Penrose inverse of the following 
two-variable polynomial matrix A(x,y): 



A(x,y) 



l-3x 5 + 9x - lOy 16 + 8x + 2y 
-7 + 9x-8y 8 + 5x~y 4 + 2x + 3y 
7-x-8y 16 - 2x - 6y -3 - 2x - Ay 



with respect to the following matrices M{x,y) and N(x,y): 
M(x,y) = 



-20-x-x 
-8- Ax- 7x 
~2{8 + Ax + 3x) 



-8-7x- Ax 
-20 + 7x + 7x 
-2(x - 5x) 



-2(8 + 3x + Ax) 

2(5x - x) 
7(-2 + x + x) 



N(x,y) 



16 + 7x + 7x 
7 — 2x — 6x 
6 — 3a: — 103; 



7 — 6x — 2x 
-2(3 + 5x + 5x) 
-2(6 + 3x + Ax) 



6 — lOx — 3x 
-2(6 + 4x + 33T) 
~3(~6 + x + x) 



The obtained weighted Moore-Penrose inverse is: 

X(x,y) = Al I(xy)N(xy) (x,y) = (60x 3 -5yx 2 -5A0x 2 + 51yx + 779x-42y-A35)~ 1 

-5x 2 + 51x - 42 -3x 2 + 8a; - 13 -3a; 2 + 33a; - 4 

-30a; 2 + 71a; + 15 42a; 2 - 5yx - 33a; + y + 15 -18a; 2 - 63a; + lOy + 105 
-2 (I0a; 2 -19a; + 12) 2 (l8a; 2 -22/a;-29a; + 2y + 17) -24a; 2 + j/a;+42a;-2j/-23 
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Let us notice that degrees of intermediate results in algorithms 4.1 and 4.2 are 
much greater than the degrees of A, M, N and X (maximum degree in this example 
are 874 and 122 of the variables x and y respectively). This is the reason why 
the algorithms for computing the weighted Moore-Penrose inverse for polynomial 
matrices are very slow (working time of the function WPolyEff for last example 
is 172.922 seconds). As we will see in the sequel, when matrices A, M and N 
are sparse, corresponding intermediate results are also sparse. Therefore, sparse 
structures introduced in the previous section improve the working time of the im- 
plementation. 

Algorithm 4.1 is tested on several random generated test examples. We tested 
variants of algorithm 4.1 using Ef and Eff sparse structures separately. In this test, 
matrices A(S), M(S) and N(S) were complex polynomial matrices of one variable 
s (i.e. holds S = (s,s)). 

We made testing for two different classes of matrices: sparse and dense. The 
measures representing sparsity of a given polynomial matrix are the same as in |12) 
(definitions 6.1 and 6.2). We are now restating these two definitions and generalizing 
them to the multi- variable complex polynomial matrices. 

Definition 5.1. For a given matrix A(S) — [ciij(S)] G C[S] rnxn (polynomial or 
constant), the first sparse number spi(A) is the ratio of the total number of 
non-zero elements and total number of elements in A(S): 

spi{Ais)) = \iMhd^M. 



The first sparse number represents the density of non-zero elements and it is 
between and 1. 

Definition 5.2. For a given polynomial matrix A(S) 6 C[S] rax ™ and S = (si,...,s p ), 
the second sparse number spa(A(5)) is the following ratio: 

, A ,^ #{(i,j,ki,...,k p ) | 0<kj<de Es A(S), Coef(a« # • • • s k /)^0} 

Sp2{A(b))= -. : -. . 

deg Si A---deg Sp ^-m-n 

By Coef (P(S), s^ 1 ■ ■ ■ Sp p ) we denoted the coefficient corresponding to s^ 1 ■ ■ ■ Sp p in 
polynomial P(S). 

The second sparse number represents density of non-zero coefficients contained 
in elements aij(S), and it is also between and 1. 

Results are presented in the next table (column d states for the degree of corre- 
sponding matrix polynomials A(S), M(S) and N(S)): 
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m 


n 


d 


Alg 4.1 
with Ef 


Alg. 4.1 
with Eff 


2 


2 


1 


0.14 


0.188 


2 


2 


2 


0.65 


1.24 


2 


2 


3 


1.92 


3.93 


3 


3 


1 


1.34 


1.32 


3 


3 


2 


9.01 


11.81 


3 


3 


3 


34.39 


48.13 


4 


4 


1 


7.87 


6.74 


4 


4 


2 


69.31 


64.48 


4 


4 


3 


461.07 


594.98 


5 


5 


1 


49.13 


58.48 


5 


5 


2 


309.38 


330.32 



s Pl {A(S)) = 0.9, sp 2 {A(S)) = 0.9 



TO 


n 


d 


Alg 4.1 
with Ef 


Alg. 4.1 
with Eff 


2 


2 


1 


0.06 


0.89 


2 


2 


2 


0.25 


0.46 


2 


2 


3 


0.60 


1.23 


3 


3 


1 


0.47 


0.68 


3 


3 


2 


4.60 


7.18 


3 


3 


3 


14.89 


24.65 


4 


4 


1 


6.10 


6.18 


4 


4 


2 


34.95 


39.68 


4 


4 


3 


256.31 


299.61 


5 


5 


1 


30.85 


39.43 


5 


5 


2 


246.32 


283.12 



a Pl (A(S)) =0.7, s P2 (A(S)) = 0.5 



TO 


n 


d 


Alg 4.1 
with Ef 


Alg. 4.1 
with Eff 


2 


2 


1 


0.04 


0.112 


2 


2 


2 


0.11 


0.263 


2 


2 


3 


0.422 


1.303 


3 


3 


1 


0.281 


0.972 


3 


3 


2 


1.367 


3.505 


3 


3 


3 


5.808 


18.449 


4 


4 


1 


1.613 


5.549 


4 


4 


2 


12.134 


27.113 


4 


4 


3 


55.139 


107.27 


5 


5 


1 


7.475 


13.582 


5 


5 


2 


84.712 


139.681 



TO 


n 


d 


Alg 4.1 
with Ef 


Alg. 4.1 
with Eff 


2 


2 


1 


0.032 


0.105 


2 


2 


2 


0.069 


0.190 


2 


2 


3 


0.187 


0.713 


3 


3 


1 


0.185 


0.675 


3 


3 


2 


0.628 


2.944 


3 


3 


3 


1.031 


3.275 


4 


4 


1 


0.987 


4.344 


4 


4 


2 


6.087 


25.263 


4 


4 


3 


27.466 


176.581 


5 


5 


1 


3.294 


15.853 


5 


5 


2 


42.159 


171.416 



s Pl (A(S)) = l, s P2 (A(S)) =0.2 



s Pl {A{S)) =0.2, s P2 (A(S)) =0.2 



All presented processor times are in seconds and the sparse numbers for matrices 
M{S) and N{S) are the same as corresponding sparse numbers for A(S). Every 
processor time is obtained by averaging working times of 15 different randomly 
generated test cases. Testing was done on Intel Pentium 4 processor at 2.6GHz 
and MATHEMATICA 5.2. We can notice that Algorithm 4.1 with an Ef structure 
showed best timings on all test cases. We have already mentioned that an Ef 
sparse structure is already implemented in MATHEMATICA. In the implementation 
we used standard built-in operators for manipulation with matrices in Ef structure. 

The first table (when s P i(A(S)) = sp 2 (A(S)) = 0.9) corresponds to dense ma- 
trices. In this case, sparse structures are not so effective because there are a lot of 
non-zero elements in all matrices and non-zero coefficients in polynomials. But we 
can notice significant improvement in working time when is applied Ef structure 
against the case when Eff structure is applied. This difference mainly comes from 
the fact that Ef structure is implemented by MATHEMATICA built-in operations. 
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The second case (when spi(A(S)) = 0.7 and sp2(A(S)) = 0.5) represents sparse 
matrices. We can notice that working times are significantly less than in the first 
case. Also here Ef structure produces less working times than Eff. 

In the third and fourth case (when spi(A(S)) = 1 and sp2(A(S)) = 0.2, and 
spi(A(S)) = sp2{A(S)) — 0.2, respectively) we deal with matrices whose entries are 
very sparse polynomials. Moreover, in the fourth case we work with matrices with 
only few non-zero elements. In the fourth case, smallest average working times are 
obtained for all considered matrix dimensions and degrees. Also we can notice that 
as sparse numbers decrease, the average working times also decrease (for constant 
matrix dimensions and degree). This holds for both sparse structures and verifies 
the theoretical results about sparse structures Ef and Eff in practice. 

We also considered simpler case: when all input matrices (A(S), M(S) and 
N(S)) and variables si, . . . , s p are assumed to be real. In that case we have only 
p variables and conjugate-transpose operation reduces only to transpose. We also 
should suppose that matrices M(S) and N(S) are symmetric in that sense. Al- 
gorithms 4.1 and 4.2 remains the same except we should change the definition 
of conjugate-transpose operator (also the implementations in MATHEMATICA). This 
case is considered in [24] and algorithms 4.1 and 4.2 are an effective versions of cor- 
responding algorithms 3.1 and 3.2 in [24]. Here working times of the algorithms are 
significantly less, and also the inverses has much smaller degrees. Results obtained 
in this special case are presented in the following table: 



m 


n 


d 


Alg 2.1 


Alg 4.1 
with Eff 


Alg. 4.1 
with Ef 


Alg 3.1 
from [24] 


3 


3 


1 


0.32 


0.23 


0.10 


0.94 


3 


3 


2 


0.69 


0.57 


0.20 


1.32 


3 


3 


3 


0.82 


1.17 


0.43 


1.84 


3 


3 


4 


1.19 


2.15 


0.73 


2.38 


4 


3 


1 


0.76 


1.26 


0.14 


1.29 


4 


3 


2 


1.29 


0.65 


0.31 


2.12 


4 


3 


3 


2.14 


1.32 


0.59 


2.42 


4 


3 


4 


2.84 


2.26 


1.01 


2.93 


5 


5 


1 


3.48 


1.45 


1.01 


3.56 


5 


5 


2 


5.90 


4.54 


2.92 


4.92 


5 


5 


3 


9.18 


8.79 


6.82 


8.27 


5 


5 


4 


12.15 


15.87 


10.85 


10.34 


6 


6 


1 


7.98 


2.65 


2.17 


8.16 


6 


6 


2 


12.93 


8.20 


7.31 


11.32 


6 


6 


3 


21.76 


18.29 


13.53 


19.42 



s Pl (A(S)) = 0.7, sp 2 (A(S)) = 0.7 



It can be seen from the table that here in all cases Ef structure was better than 
Eff (both with using Algorithm 4.1). Both effective algorithms was significantly 
better than Algorithm 2.1 (for rational matrices) and Algorithm 3.1 from [24]. For 
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smaller values of d, Algorithm 2.1 was better than Algorithm 3.1 from [53] due to 
the implementation details. 

All presented results leads us to the same conclusion: the best choice for com- 
puting weighted Moore-Penrose inverse for polynomial matrices is Algorithm 4.1 
with the sparse structure Ef . 

6 Conclusion 

We extend the algorithm for computing the weighted Moore-Penrose from [20] to 
the set of multiple-variable rational matrices with complex coefficients. We adapt 
previous algorithm to the set of polynomial matrices. We consider two effective 
structures which make use of only nonzero addends in polynomial matrices and im- 
prove previous results on the set of sparse matrices. In the last section we presented 
an illustrative example and compared various algorithms. 
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